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Abstract 

The most widely-used density functionals for the exchange-correlation energy are inexact 
for one-electron systems. Their self-interaction errors can be severe in some applications. The 
problem is not only to correct the self-interaction error, but to do so in a way that will not 
violate size-consistency and will not go outside the standard Kohn-Sham density functional 
theory. The solution via the optimized effective potential (OEP) method will be discussed, first 
for the Perdew-Zunger self-interaction correction (whose performance for molecules is briefly 
summarized) and then for the more modern self-interaction corrections based upon unitarily- 
invariant indicators of iso-orbital regions. For the latter approaches, the OEP construction is 
greatly simplified. The kinetic-energy-based iso-orbital indicator r ( ^(r)/r cr (r) will be discussed 
and plotted, along with an alternative exchange-based indicator. 

1 Introduction 

Paradoxically, some of the most widely-used and reliable theoretical approaches to many-electron 
systems, including the local spin density [Tj (LSD) and generalized gradient 2 (GGA) approxi- 
mations, are not exact for one-electron systems. The error they make in these systems is called 
the self-interaction error. An early self-interaction correction [3J @] (SIC) to LSD led to dramatic 
successes and failures, and has been largely bypassed by the development of GGA, which provides a 
more uniform improvement over LSD and has an easier Kohn-Sham theory implementation. Since 
the time finally seems ripe for the development of reliable self-interaction-free approximations, this 
article will discuss some of the possibilities for self-interaction correction within Kohn-Sham density 
functional theory pQ. 
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In Kohn-Sham theory, the many-electron ground-state spin densities Wt(r) and n\ (r) and energy 
E are predicted by self-consistent solution of the one-electron equations 



(p aa (r) = £ acT ip a a(r), (1) 



2m 

n a (r) =^2\p aa \ 2 Q(fi- e aa ) (2) 

a 

where the spin-dependent effective potential is 

<(r) = v a (r) +e 2 J d\> 0^ + < c (r). (3) 

In Eq. ©, v a (r) is the external potential created by the nuclei and external scalar fields, n = n-^+ni 
is the total electron density, and 

5E xc [n hni ] 

xc{ ' 5n a (r) [ > 

is the exchange-correlation potential, is the Fermi level, and Q(x) is the step function (0 = 1 for 
x > and = for x < 0.) The energy is 

E = T s [n h n i \+Y J J d 3 rn a (r)v a (r) + j J d 3 r J d 3 r' n ^^ + E xc [n h n;]. (5) 

The non-interacting kinetic energy is 

T s [n T ,nJ =2 / ^s(r) = £ / ^(r), (6) 



where 

k(r) = ( "^ y2 ) faa(r)@^-s aa ), (7) 



fi 2 



ft 2 



^(r) = E— |V^ CT (r)| 2 0(/i- £c 
2m 



a 

,2 



= t CT (r) + ^V 2 n (7 (r). (8) 

Since the Kohn-Sham orbitals vw( r ) are functionals of the spin densities rtt and jii , so is T s . 

While the equations of the previous paragraph are exact in principle, in practice E xc [n^,n^] = 
_E x [re|,nj] + 2£ c [n-f ,njj must be approximated. In the local spin density approximation pQ, 

E^ D [n hni ] = J d 3 rnelf(n hni ), (9) 

where e™ is the known exchange-correlation energy of an electron gas with uniform spin densities 
n-f, n[. In the generalized gradient approximation [2], 

,-iGGAu ™ 1 _ /" j3„_GGA 



C A K.»J = / Aner^t^l.Vnt.Vnj). (10) 
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These approximations are exact for a uniform density and accurate for a slowly- varying n(r), but 
are not exact for one electron densities. 

One-electron densities are fully spin-polarized (e.g., nt = n and ny = 0). In a one electron 
system, the exchange energy must cancel the Hartree electrostatic energy: 

E xM = -^[^r^^ (N = N, = 1) (11) 
2 J |r — r'| 

and the correlation energy must vanish: 

E c [n h 0}=0 (jV = iV T = l). (12) 
To satisfy Eqs. (|llj) and (|12l) . Perdew and Zunger [5] proposed a self-interaction correction to LSD: 

^ = ^ SD h,n x ]-g^ (13) 
^ C ' 0<T W = ^ C SD ' CT (K, «J; r) - e 2 / rfV - ^([n^, 0]; r), (14) 

where 

n QCT (r) = |v3j CT (r)| 2 9(/i - e aa ) (15) 
is an orbital density. The potential (|14j) has the correct asymptotic behavior 

e 2 

v xc( r ) ~> as r —> 00 (16) 

r 

as one moves away from any compact system, while u^ D,CT (r) tends to zero exponentially in this 
limit. Unlike LSD, the SIC exchange-correlation energy displays 0] a derivative discontinuity very 
much like that of the exact E xc [nun\]. 

There is no unique way to make a self-interaction correction, and alternatives to Eqs. (|13|) - 
(|15|) have been proposed [3 El HI El E] • But Eqs. (|13|) - Q15|) have been widely tested for atoms 
UU1 EU d ECU , atomic ions UH , molecules [13111117111311111101111111111121121113 and solids 
|201 I27j (see earlier references in Ref. SIC is exact for one-electron systems, and usually 

accurate for strongly localized electrons. For covalent molecules near equilibrium, it has been 
argued [13 H3 H3 that the self-interaction error in LSD and GGA exchange f,28 4 is needed to 
mimic the effect of static correlation on the electron density. 

In applications to molecules, the performance of SIC is somewhat mixed. Total energies are 
better than in LSD, and the highest occupied orbital energy is much closer to minus the ionization 
potential than in LSD jlbl I18j . The localized SIC valence orbitals correspond to the localized 
bonds and lone pairs of chemical intuition ^3- SIC significantly improves the energy barriers to 
chemical reactions |17l I22j . but net reaction energies are less strongly improved relative to LSD 
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|22| . Many nuclear magnetic resonance properties of molecules are improved by SIC |21j . There 
are relatively few studies of atomization energies in SIC, but there seems to be an improvement 
over LSD for the cases studied: LI2 ^H], O2 (UJ, and N2 PH. The most disappointing results are 
the SIC bondlengths, which are shorter than the experimental ones by ~ 0.07 bohr on average |18j . 
while the LSD bond lengths are much more realistic. 

2 Perdew-Zunger SIC within Kohn-Sham Theory 

The SIC of Eqs. ()13|) - Q15|) goes outside the Kohn-Sham scheme by introducing an orbital- dependent 
effective potential vf a (r). As a result, the self-consistent SIC orbitals are not Kohn-Sham orbitals, 
and are not even strictly orthogonal unless off-diagonal Lagrange multipliers are introduced. The 
SIC orbitals tend to localize around atomic centers, while the Kohn-Sham orbitals are delocalized 
canonical or molecular orbitals. The SIC orbitals can be found, even for molecules, by directly 
minimizing Eq. (|13j) under the constraint of orbital orthogonality |151 1161 1231 1241 125j . 

Although not a Kohn-Sham theory, the Perdew-Zunger SIC belongs [3] to a wider class of density 
functional theories. At least to the extent that the SIC orbitals are localized, it is also a size- 
consistent theory jl], i.e., one which works consistently well for small or large systems. 

But there are clearly computational and conceptual advantages to Kohn-Sham theory, not only 
for the ground state but also for time-dependent processes and excitations. To bring SIC under 
the umbrella of Kohn-Sham theory, one must construct a common effective potential for all the 
occupied orbitals of spin a. Especially in the context of time-dependent DFT, different procedures 
emphasizing computational simplicity have been suggested to construct a common local potential 
[213 EUJ , and the influence of the self-interaction correction on optical properties of atoms |3T?j . 
the molecule N2 [H2| and clusters jlHl BH] has been discussed. A rigorous way of constructing a 
common potential is given by the optimized effective potential (OEP) method |341 1351 137)] . For 
any orbital functional £[{i^ aff }], the OEP method delivers a Kohn-Sham potential and a set of 
Kohn-Sham orbitals which minimize that functional. When the orbital functional is Hartree-Fock, 
there is no problem, but when it is SIC (Eq. (|13j)) the resulting scheme is not size-consistent: 
Applied to one atom, where all the Kohn-Sham orbitals are localized, this scheme will deliver a 
properly self-interaction-corrected energy. But, applied to a periodic lattice of atoms separated by 
large lattice constants, where all the Kohn-Sham orbitals are delocalized, this scheme will produce 
no self-interaction correction to the energy of an atom, since the sum in Eq. (|13j) will then vanish 
on a per-atom basis [3]. The considerations put forward in Refs. 013] suggest that this is true 
for all Perdew-Zunger-like SIC schemes that directly use the Kohn-Sham orbitals. Therefore, such 
schemes would be good for atoms, but would degrade for molecules or clusters as the number of 
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atoms increased. 

A clue to the solution of this problem was given in the work of Pederson, Heaton and Lin |16| . 
who introduced two sets of occupied orthonormal orbitals related by unitary transformation: the 
localized SIC orbitals, and the delocalized canonical orbitals. Garza, Nichols and Dixon j2U] pro- 
posed that the canonical orbitals could be Kohn-Sham orbitals belonging to an optimized effective 
potential Vg(r) constructed from the localized orbitals. In their work, and in that of Patchovski 
and Ziegler the Krieger-Li-Iafrate approximation [37j to OEP is used, as is a standard 

(non-optimal) localizing transformation. 

As an exactification of this approach, the correct Kohn-Sham version of Perdew-Zunger SIC 
would be conceptually this: Start with a given external potential v a (r) and electron number N. 
Form a trial effective potential uj(r), and solve Eq. (^Q) to find the corresponding occupied Kohn- 
Sham orbitals. Then find the unitary transformation to localized orbitals that minimizes Eq. (|13j) . 
Finally, choose the effective potential that delivers the lowest minimum of Eq. (|13l) , 

3 Unitarily Invariant Iso-Orbital Indicators 

The prescription outlined above for the implementation of the Perdew-Zunger self-interaction cor- 
rection to LSD (or GGA) within Kohn-Sham theory was greatly complicated by the fact that the 
self- interaction correction was not invariant under a unitary transformation of the occupied orbitals. 
This section will discuss self- interaction corrections that are unitarily invariant, and thus can be 
implemented within Kohn-Sham theory by a direct application of the OEP method |34l 1351 136| to 
the Kohn-Sham orbitals. This subject is timely because of the recent appearance of accurate and 
efficient solutions EH] to the OEP problem. 

A Slater determinant of occupied orbitals of a given spin a is invariant under unitary transforma- 
tion of those orbitals, and so is any quantity that can be constructed from the Slater determinant, 
such as the spin density of Eq. @ or the kinetic energy densities of Eqs. Q) and (|HJ). The one- 
electron density matrix 



is also invariant. (The step function must of course be re-interpreted as a restriction to the occupied 
orbital space.) The exact exchange energy 



p a (r,r 7 ) =^2(p* a(T (r)<p aa (r')@(fi - s aa ) 



(17) 



a 




(18) 



is clearly invariant, as is the the local exchange energy per electron e x (r): 




(19) 
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On the "Jacob's Ladder" |38| of density functional approximations, full freedom from self- 
interaction error is achieved only at the hyper-GGA level, which employs full exact exchange and 
a highly nonlocal functional of the occupied orbitals for correlation. A somewhat different way to 
eliminate the self- interaction error is via a local hybrid functional 39, . But in either case one needs 
an iso- orbital indicator to identify regions of space in which the electron density is dominated by a 
single orbital shape. The iso-orbital regions where n^n^ = are one-electron regions in which the 
correlation energy per electron e c (r) can and should be zeroed out by a self-correlation- free density 
functional. 

The exact exchange energy of Eq. (|18|) is self-interaction free, since for a one-electron (N = N a = 
1) ground-state p CT (r, r') = n l J 2 (y) n l J 2 {y'). Thus 

1/2, \ 1/2, a 

is an iso-orbital indicator which equals unity when both r and r' are in an iso-orbital region. 
However, as r' — > r, y a (r, r') tends to 1 in any region, iso-orbital or not. This problem does not 
arise for 

= lim = ^W. (21) 

V ' r'-r V r - V r ,/v(r,r') r ff (r) V ; 

Eq. (|21|) provides a point-by-point iso-orbital indicator which equals unity in any iso-orbital region 
and is otherwise bounded between and 1 jlU]. In Eq. 1)211) . T a (r) is the kinetic energy density of 
Eq. (JHJ), and 

T „ _ pEagf (22) 

8m n a {r) 

is the von Weizsacker or bosonic kinetic energy density. For a uniform density, x a (r) vanishes 
everywhere. 

x a (r) of Eq. (|21|) is clearly invariant under unitary transformation of the occupied orbitals. The 
idea of using the condition T CT (r) = r ( ^ / (r) to identify an iso-orbital region and zero out the self- 
correlation goes back to Colle and Salvetti [U, but in density functional theory to Becke [HI an d 
Dobson 43 . x CT (r) is an ingredient of self-correlation free meta-GGA's including those of Refs. 
and and of local hybrids 39 and hyper-GGA's [SHI- 

Fig. ^ shows x,j(r) of Eq. (|21j) . plotted as a function of the distance from the center of two 
spherical densities, the Ne atom and a jellium cluster with 34 electrons having the average valence 
electron density of Na (r s = [3/(47rn)] 1 / 3 = 3.93 bohr). ^(r) correctly identifies the density tails 
as iso-orbital. In the Ne atom, the Is core is also found to be nearly iso-orbital. The interior of 
the jellium cluster, however, is found to be a region of strong orbital overlap, as expected. The 
densities and orbitals have been evaluated by solving [^31 EH] the OEP problem for exact exchange. 
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Figure 1: The kinetic-energy-based iso-orbital indicator x(r) = t w (r)/r(r) of Eq. (|21[) . as a function 
of the radial coordinate r (in bohr) for two spherical densities: the neon atom, and a jellium cluster 
resembling Na34, with a radius of 12.7 bohr. 



While x a {r) of Eq. (j21|) seems to be a satisfactory iso-orbital indicator, it does display an order- 
of-limits problem Define a = (r — t w )/t , where To = ^;(3vr 2 ) 2 / 3 n 5//3 , and p = | (t w /t) = 
|Vn| 2 /[4(37r 2 ) 2 / 3 n 5 / 3 ]. Then, for n T = n b 



a; 



lim lim x 



1 + i-' 

5 p 



1 but 



lim lim x 



0. 



(23) 
(24) 



This problem shows up in nearly-iso-orbital (a — > 0) regions where the gradient of the density 
approaches zero (p — > 0), and thus perhaps at covalent bond centers. 

Because of the order-of-limits problem of t w /t, it may be worthwhile to consider alternative 
iso-orbital indicators. For example, the exact exchange potential v x (r) and the exact exchange 
energy per electron e x (r) (Eq. I[19|)) of a spin-unpolarized system are related in the iso-orbital limit 



N 


ratio 


N 


ratio 


2 


1.000 


40 


0.547 


8 


0.834 


58 


0.528 


18 


0.679 


92 


0.528 


20 


0.640 


106 


0.390 


34 


0.571 


138 


0.442 



Table 1: The ratio of Eq. (|28|) for closed-shell jellium spheres of increasing electron number N. 
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Figure 2: The exchange-based iso-orbital indicator x 
densities as in Fig. ((J). 



Vex 



4) of Eq. (|27j). for the same 



by 



and in the uniform-density limit by 

v x (r) 



e x (r) 

4 
3 



(TV = 2), 



(uniform density). 



One might define 



x(r) 



^x(r) 
e x (r) 



(25) 



(26) 



(27) 



as an alternative iso-orbital indicator, which varies from 1 in the iso-orbital limit to in the uniform 
limit. Tabled shows that 

(28) 



J d 3 rn(r)v x (i 



J d 3 rn(r)e x (r) 3 

varies almost smoothly from 1 for the N = 2 jellium cluster to 0.4 for the largest cluster studied 
here. Fig. [21 however shows that x(r) of Eq. (|27|) can be negative, fails to recognize the Is core of 
the Ne atom as a strongly iso-orbital region, and fails to recognize the interior of the jellium cluster 
as a region of strongly overlapped orbitals. 

Thus it seems that x a (r) of Eq. (|21jl is the preferred iso-orbital indicator. Note that Eq. (|22jl 
assumes that the orbitals can be chosen to be real. If the state of interest has a nonzero current 
density 



7; 



occup. 



j CT (r)=Re— J2 ¥C( r )-Vv?a<x(r) 



then H3 HH 



rf(r) 



h 2 |V?y(r) 
8m no-(r) 



+ 



m\j a (r)\'- 
2n ff (r) 



(29) 



(30) 
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To ensure that x a equals one for a one-electron density and zero for a uniform density (with or 
without a uniform current), the best choice may be x a = jf a where = — m|j t7 | 2 /(2n tT ) 
and f a = r a — m|j cr | 2 /(2n (7 ). In this way, the self-correlation error can be corrected even in a general 
excited state. 

4 Conclusions 

For many standard applications of ground-state density functional theory, the self-interaction er- 
rors of modern GGA's and meta-GGA's are relatively benign. There are a few striking exceptions 
to this rule, such as the binding properties of diatomic molecules with an odd number of valence 
electrons |191 I47j and the static (hyper-) polarizabilities of long-chain molecules |48| . For appli- 
cations involving time-dependent and excited-state Kohn-Sham density functional theory |32l I49j . 
the self-interaction errors can be severe. 

While the Perdew-Zunger self-interaction correction to the local spin density approximation can 
now be brought under the umbrella of Kohn-Sham theory, the development of more sophisticated 
functionals and optimized effective potential methods suggests that general-purpose self-interaction- 
free density functionals will be developed soon and implemented within Kohn-Sham theory. Such 
functionals may well include full exact exchange plus highly nonlocal correlation based in part upon 
unitarily- invariant iso-orbital indicators such as r < ^(r)/r cr (r). 

Acknowledgements: Our investigation of Eq. (|27|1 was triggered by discussions with Prof. Dietmar 
Kolb. S.K. acknowledges financial support by the Deutsche Forschungsgemeinschaft under an 
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